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ABSTRACT 

Context. An analytical solution of the generalized diffusive and convective transport equation is 
derived to explain the transport of cosmic ray protons within elliptical galaxies. 
Aims. Cosmic ray transport within elliptical galaxies is an interesting element in understanding 
the origin of high energetic particles measured on Earth. As probable sources of those high ener- 
getic particles, elliptical galaxies show a dense interstellar medium as a consequence of activity 
in the galactic nucleus or merging events between galaxies. Thus it is necessary for an appropri- 
ate description of cosmic ray transport to take the diffusive and convective processes in a dense 
interstellar envirorraient into account. Here we show that the transport equations can be solved 
analytically with respect to the given geometry and boundary conditions in position space, as well 
as in momentum space. 

Methods. From the relativistic Vlasov equation, which is the most fundamental equation for a 
kinetic description of charged particles within the interstellar medium in galaxies, one finds a 
generalized diffusion-convection equation in quasilinear theory. This has the form of a 'leaky 
box' equation, meaning particles are able to escape the confinement region by diffusing out of 
the galaxy. We apply here the 'diffusion approximation', meaning that diffusion in gyrophase and 
pitch angle are the fastest particle-wave interaction processes. An analytical solution can be ob- 
tained using the 'scattering time method', i.e. separation of the spatial and momentum problems. 
Results. The spatial solution is shown using a generalized source of cosmic rays. Additionally, 
the special case of a jet-like source is illustrated. We present the solution in momentum space with 
respect to an escape term for cosmic ray protons depending on the spatial shape of the galaxy. For 
a delta-shape injection function, the momentum solution is obtained analytically. We find that the 
spectral index measured on Earth can be obtained by appropriately choosing of the strength of 
Fermi I and Fermi II processes. From these results we calculate the gamma-ray flux from pion 
decay due to proton-proton interaction to give connection to observations. Additionally we deter- 
mine the escape-spectrum of cosmic rays. The results show that both spectra are harder than the 
intrinsic power-law spectrum for cosmic rays in elliptical galaxies. 

Key words, cosmic rays: general - transport model ~ diffusive processes - convective processes 
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1. Introduction 



Since their discovery by Viktor Hess in 1912, cosmic rays have been one of the biggest fields of 
interest in astrophysics, and yet the origin of these particles is still an open question. Fully ionized 
atomic nuclei reach the Earth coming from outside the solar system with very high energies up 
to 10^*^ eV. Most of them with energies < 10'^ eV seem to originate in the Milky Wa y, while the 
highe st energetic ones are considered to have an extragalactic origin (for a review see 



Hoerandel 



2007) 



The most accepted model for the origin of ultrahigh ene rgy cosmic rays (UHECRs) is accel- 
eration in shock fronts due to Fermi processes (lFermilll949h : hence , the main source s are gamma 
ray bursts (GRBs), active galaxies like active galactic nuclei (AGN) ( Tavecchiol2005 ). or colliding 
galaxies. The last two are specially interesting for two reasons. First UHECRs can be generated 
in elliptical galaxies. Second the increase in the interstellar me dium density in objects of these 
types influences cosmic ray transport (see lBekki & Shiovalll998l) . But since all elliptical galaxies 



have an interstellar mediu m due to star winds, we conclude that the study of transport processes is 
interesting in general (see Knapp|[l999 ). 

In particular, as a res ult of the GZK-e ffect, very close AGN are the most probable candidates 
for UHECR sources (see Biermannlll995 ). One of these nearby AGN is the giant elliptical galaxy 



M87. It is propose d that this galaxy is respons i ble for acceleration of cosmic ray protons due to 



Fermi I processes ( Blandford & Ostriker 



1978 



dReimer et al 



Rieger et al J 120071) in s hockfronts within the jet 



20041) . For an overview of proton acceleration in jets (see 



Mannheim 



1993). Since 



particle acceleration in a jet is located within active galaxies surrounded by an interstellar medium, 
the high energetic protons undergo physical transport processes before they escape out of the galaxy 
and reach the detectors on Earth, making a model for cosmic ray transport in elliptical galaxies 
inevitable. This helps give an answer to the question about the origin of cosmic rays. 

Progress has been made in the field of modeling cosmi c ra y transport with the numerica l 
description of transport processes bv lOwens & Jokipiil (119771) and IStrong & Moskalenko (1998). 
Never t heless we follow the ba s ic idea s pre sented in the underlying papers of lLerche & Schlickeisei 
J 19851) . Iwang & SchUckeisej Jl987n . and Lerche & Schlickeisej Jl988l) . who used an analytical 
description of cosmic ray transport. In relation to our work such, a treatment has the following 
advantages: the model is adequate for cosmic ray transport within any kind of elliptical galaxy 
including arbitrary cosmic ray sources and the physical parameters involved in our model can be 
easily fitted to measurements. After all, our analytical model can serve as a test case for more 
profound numerical models. 

In this paper we solve the cosmic ray transport equation analytically with respect to a kinetic 
description of the interstellar plasma in elliptical galaxies. Special attention is paid to the spatial 
transport of charged nuclei. In addition, the solution of the momentum equation is derived to ex- 
plain general properties of this model. As a result, we present illustrative examples of spatial, as 
well as momentum, cosmic ray transport for given sources of charged nuclei. To show the connec- 
tion to observations, we calculate the gamma-ray flux from neutral pion decay. These mesons are 
produced by inelastic scattering processes between cosmic ray protons. The resulting power-law 
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spectrum is slightly harder than the intrinsic one for cosmic rays. Finally, we present the escape 
spectrum of charged particles leaving elliptical galaxies. Similar to the gamma-ray flux, this spec- 
trum is flatter than the intrinsic one. 



2. Basic equations 



To describe the propagation of cosmic ray nuc l ei wi thin elliptical galaxies, we follow 



Lerche & SchUckeised d 19851) 



Lerche & Schlickeisen (119881) and 



Schlickeiseij (120021) . For the de- 



scription of transport processes they use the 'diffusion approximation', which means that the fastest 
particle-plasma wave interaction processes are d iffusion in gyrophase and p itch angle. Thus fol- 



lowing |Joki2ii| (Il96q) 



Hasselmann & Wibberena (119681) . and 



Skillina (119751) . we take an isotrope 



particle distribution function in momentum space. Here we ideaUse the interstellar medium as a ho- 
mogeneous volume containing primary cosmic rays being accelerated from the thermal background 
medium and secondaries result ing from fragmentation of prim aries having a negligible abundance 
in the background medium (cf. 



Havakawalll969: 



I 



Cowsiklll98C). 



The transport of these particles at lar ge momenta (p > 10 GeV c ' nucleus ') is described 



Schlickeise Such a treatment is suitable to short 



by the steady-state transport equation (e.g. 
timescales of diff'usive and convective processes compared to the dynamical timescale of the galaxy 
(fdyn ~ 10^ years). This is true in the case of high energetic particles. We assume a spatial diffusion 
coefficient K(r) of 10^"^ cm^ s ' at p = 1 GeV, which is slightly larger than the value measured in 



the Milky Way {K(r) = 10" cm , (cf . ISchlickeiseiil2002h because of diffusive processes being less 



effective in elliptical galaxies. Consequently we get for protons with TeV-energy a timescale of 
X 10^ years. Furthermore, the dynamical age of the galaxy has to be greater than the timescale of 
source variability to obtain an appropriate description. This is usually given, since the size of the 
accretion region onto the central black hole is of the order of a few light-days so that a maximal 
variability timescale of some days is assumed. 

At large momenta, spatial diffusion in turbulent magnetic fields dominates convection in the 
galactic wind, so that we find a transport equation for the phase space density /(r, p) in spatial 
coordinates r and in the momentum coordinate p: 



J:rf + -Cpf + S(r,p)^0. 



The spatial operator Xr is defined by 



£Ar,p) = V[K(r,pW] 



(1) 



(2) 



containing spatial diffusion with the spatial diffusion coefficient K{r, p) - K{r)K{p), where K{p) 
denotes the dimensionless dependence on the momentum variable p without loss of generality. The 
momentum operator 



A(r,p).p-^^ 



1 



(3) 



describes momentum diffusion by second-order Fermi processes (D(p)), energy gain due to first 
order Fermi processes (pgain). as well as conti r iuous (pioss) and catastrophic (tc) momentum loss 
processes. As shown in lLerche & Schlickeisen (119851) . fully-ionized particles heavier than protons 
have the same Fermi acceleration rates as protons, so hereinafter momentum means momentum per 
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nucleon. We are interested primarily in the behaviour of the cosmic ray primary spectrum, so we 
do not take secondary particles due to fragmentation of primaries into account. On the other hand, 
we allow fragmentation of primaries as a general loss process. 

A link between spatial and momentum diffusion processes can be seen in the relation between the 

two diffusion coefficients ^ 

D(r,p)^ ^'^ , (4) 
" K(r,py 

where C\ stands for the proportionality factor being independent of r and p. This close connection 
arises from the same basic physical process behind spatial and momentum diffusion: Protons are 
scattered in pitch angle due to the magnetic fields of MHD plasma waves causing spatial diffusion 
along ordered magnetic field lines, whereas cyclotron damping of the electric field associated with 
MHD waves affects diffusion in momentum space. For that reason it is necessary to solve this 
model in spatial coordinates as well as in momentum coordinates to get an acceptable description 
of transport processes in elliptical galaxies. 



3. 'Scattering time' method 



We use the 'scattering time' method proposed bv lSunyaev & Titarchukl(ll98C ) to g et an important 



class of exact analytical solutions of Eq.([T]) following 
that the spatial and momentum operators can be separated as 

£r(r,p)^h(p)Or(r) , £p(r, p) ^ g(r)Op(p). 



Wang & Schlickeiseij (Il987l) . This implies. 



(5) 



For ease of exposition, the source function S (r, p) is also a product of two separable functions, i. 
e., 

S(r,p)^q(r)Q(p). (6) 

As an aside we note that the requirement Eq.(|5]) is trivially fulfilled, if K(r) is constant (^r(r) = Kq) 
and /igain, pioss, as well as r^,, are all independent of spatial variables. Thus we use the following 
model containing constant factors a„ with n = 1,2,3,4 for the proportionality factors independent 
of spatial and momentum coordinates. First we take D(p) = a2p^/K(p) to describe the momentum 
diffusion coefficient solely in momentum space. Furthermore we assume 



Pgain = aip/K{p) 



(7) 



telling that first-order Fermi acceleration is related to the (momentum dependend) spatial diffusion 
coefficient due to MHD plasma-wave-scattering interactions within the acceleration process. The 
continuous loss term /)ioss is independent of spatial coordinates leading to 



Similarly, we set 



Piossip) = -a^pip). 



Tcip) = fl4 e(p) 



(8) 



(9) 



for the catastrophic loss time due to fragmentation. 

Under these conditions, in addition to Eq.© and (|6]l, we can find the formal mathematical solution 
of Eq.([T]i as a convolution of the spatial and momentum solution functions T(r) and M{p) following 



de Freitas Pachecol (119711) : 



f(r,p) 



du T(r, u)M(p, u). 



(10) 
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where T(r, u) has to satisfy the given spatial boundary conditions, and 

dT 1 ^ 
— = —-OrT 

du g{r) 



with 



and the conditions 



a = V [KoV] = KqA , g(r) = 1 

T(r,u = 0)^q(r)/g(r)^q(r) 
T(r, M = oo) = 0. 

Here, M(p, u) has to satisfy the given spatial boundary conditions, and 

dM 1 



with 



1 d 
p^dp 



-OpM 
ou h{p) 



+ a^p pip) 



k(p) dp k(p) 



fl4 



; h(p) = k(p) 



and the conditions 



M(/7,M = 0) = Q(p)/h(p) = Q(p)/K(p) 
M(p, u - oo) -0. 



(11) 

(12) 

(13) 
(14) 

(15) 
(16) 

(17) 
(18) 



As a consequence of the formal mathematical solution, we note that we have to solve two partial 
differential Eqs.lfTTTi and (fTsT i instead of the much more complicated differential Eq.([Tli. 



4. Results 



4.1. Spatial solution 

The most convenient way to find the formal solution of Eqs. lfTTl i and (fTSI l is to start with the spatial 
problem. As can be s een from Eq.(fT2li. the spatial operator Or is of Sturm-Liouville type (cf. 
Arfken & Weben2005h and therefore has a complete eigenfunction system Ei(r). As a consequence, 
the solution function T(r, u) can be expanded in this orthonormal system as 



Here the A,(r) are defined by 



r(r,M) = ^A,(r)e-'^'". 



A,(r) = UiEiir), 



(19) 



(20) 



implying that the coefficients a, weight each eigenfunction Eiir). The /l, denote the eigenvalues of 
Or implying the special spatial geometry. 

The shape of elliptical galaxies is adjuste d to the cosmic ray transport Eq.(fT2b using prolate 
spheroidal coordinates as they are defined bv lAbramowitz & StegunI ( 119721) : 



r\ + r2 



(21) 



2/ 2/ 

As can be seen from Fig.® ri and r2 are the distances to the foci of the confocal ellipse, where 2/ 
denotes the distance between the two foci F\ and F2. Additionally, we use the variable 4> for the 
usual azimuthal dependence like in spherical coordinates. The following relations give the relation 
between these coordinates and the semi-major axis a and the semi-minor axis b, respectively: 



(22) 



6 



/ ea=f 


a \ 






\ ^^^^ 





Fig. 1. Schematical view of an ellipse with its fundamental properties. 



The numerical excentricity e has a direct relationship to the coordinate ^ via 

e = //a = l/f. (23) 

The variable r\ is defined as rj - cos 6 with 9 the angle between the line on which the foci lie, and 
an arbitrary point on the ellipse. As a result the variables are defined in the range 

^ e [l;co[ , T] e [-1;1] , e [0;2n]. (24) 




Fig. 2. Left: Gray-shaded plane cuts ellipse leading to a cut view like the graph on the left-hand 
side. 

Right : Direction of uni t vectors of the variables ^, rj, and in prolate spheroidal coordinates. 
From lweissteiii ( 1999 ). 



Figure (|2]l shows an illustration of the definition of the three spatial variables ^, rj, and 
Because of these definitions, we can write Eq. (fT2l) as 



^^0 



dT 



d_ 

di] 



J dT 

d-T?')^ 

or] 



(25) 
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Here we used the Laplacian in prolate spheroidal coordinates. The general solution can be obtained 
by consecutive separation of variables (see Appendix A): 

T{^, I], <p,u) ^ ^) ^ '?) ^ cosimcp) X exp(-^2^) = 

(2'« + '-)! ™„ 



r=0,l 



-1 n 
I ^2 _ Y 



X 



2 



X 

)■=(), 1 



2c^ 

x^?;:'"(c)P™+,(?7) X cos(m0) X exp(-A:2M)} . (26) 

Here the P"^^^j.{rj) denotes associated Legendre functions of the first kind, order m + r, and the 
-/„+ i (c^) denote Bessel functions of the first kind and of order n + j. The sum is extended over even 
values of r as indicated by the mark '. The factor c is given by c = 

To define reasonable boundary conditions, we assume a 'leaky box' model. Cosmic ray particles 
are trapped by disordered magnetic fields within the confinement region of an elliptical galaxy. In 
this they undergo diffusive and convective movements. At the edge of the box, leakage out of the 
confinement area is possible. 

As an illustrative example for spatial boundary conditions, we show the solution depending 
on a constant source function over the elliptical galaxy in Appendix B. To be more specific, we 
take a jet-like source function here. The jet points in the direction /yinj (represented by a Dirac 
delta function) with a length scale chosen to be /max for any choice of ^ being smaller than an 
arbitrary maximum value of the confinement region ^c- Particles leak out at the edge of this region. 
Such a boundary condition is known in the literature as a 'free-escape' condition. For a realistic 
assumption we decide to let the jet end smoothly (see the 'Fermi' function in Eq.(|29l)). We neglect 
any dependence on for an adequate illustration. These conditions are taken into account by 

r(^ = ^„77,0,M) = O, (27) 

by a periodical boundary condition in 77 

r(^,7?-i,0,M) = r(^,7? = +i,0,M), (28) 

and by 

^ [exp(4(/- /„,,)+ 1] 

Finally after some calculations in which Eq.(l26]l has to match the boundary conditions Eq. (l27H29] l. 
we derive the general solution (cf. Appendix B): 



00 00' ^ I £\ 



1=1 r=0 

The weighting factors are 



- . . u 



(30) 



X -T. . 



exp[4(/-/_)] + l f [y,.,3(y,>)]' 



Here the y,> stands for zeros of (;c). The sum over r is extended over even values of this param- 
eter. We see that the general solution of Eq.(|25]) with respect to spatial boundary conditions indeed 
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Fig. 3. Graphic demonstration of the cosmic ray particle density given by the solution Eq.(l30]l with 
the weighting factors OTT l. The cosmic ray particle density is normalised and given in arbitrary 
units. As an effect of chosen coordinates only 'one half of the galaxy is visible. For illustration we 
plotted T{f, Ti,u - 0) with a constant ^ = 1 .2 specifying an E4 elliptical galaxy. 

has the form of the eigenfunction expansion Eq.(fT9l) with Eq.(l20l). This solution is shown in Fig. 
([3]) with respect to 100 spatial eigenfunctions (imax - 10, r^ax = 18). The jet is responsible for the 
cosmic ray particles distributed over the whole galaxy due to diffusive processes. It is broadened 
with a larger distance to the centre. This effect is caused by the chosen geometry, whereas the loss 
of magnetic collimation in real astrophysical sources can provide this. To illustrate the shape of an 
elliptical galaxy we plotted T{ f, 77, m = 0) with a constant ^=1.2 specifying an E4- galaxy. Using 
a more complicated source function, a better physical description of particle distribution within 
elliptical galaxies can be obtained. 

4.2. Consistency checks of the formal spatial solution 

To get a better understanding of our model, we prove the formal mathematical solution (Eq.(l26]l) 
of the spatial cosmic ray transport equation. Our discussion is related to the solution found in the 
illustrative example in Appendix B, but can similarly done with Eq.(l26]i: 




(32) 



with the weighting factors 



(33) 



Here we used periodical boundary conditions for 77 and <p 



T(£,i]^ -l,4>,u) = T(^,T]^ +1,0, w), 



(34) 



T(^, T],(f>^ In, u) = r(^, /?, <^ = 0, m). 



(35) 



and 



r(^,77,0,M = O) = ^o0(^,-^). 



(36) 
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First, this solution has to accomphsh the given spatial boundary conditions. As suggested in Eq.(l36]l 
we used a constant source function over the whole size of the galaxy in order to lose all angular 
dependencies. This behaviour is caused by the source function not having any dependence in rj 
and <p, but also due to the periodic boundary conditions chosen for these variables. Such an ef- 
fect can also be seen in a spherical or disc geometry. To the spherical geometry as performed by 
(119871) . we added angular dependencies and a constant source function over the 



Schlickeiser et al 



galaxy as an example. The weighting factors of the general solution functions turn out to disappear 
in any case except for the angular separation constants being zero (the angular solution functions 
are spherical harmonics in this case), which means that there is no dependence on these variables 
as it is expected. 




Fig. 4. Normalised weighting factors a, dependent on the number of zeros y, in Eq.(l32]i 



Second, we show that the sum of solution functions (cf. Eq.(l26]l). together with the weighting 
factors in addition to the given spatial boundary conditions converges. For the illustrative example 
in Appendix B, the normalised expansion coefficients a,- are shown in Fig.©. To obtain applicable 
results it is essential to include only the first few ones depending on the requested accuracy. The 
resulting solution function is shown in Fig.©. In the special case of a constant source function, 
the solution reproduces the boundary condition 6{^g - ^) for the variable ^. We used - A for 
illustration. If we take many eigenfunctions into account, we see at the discontinuity points ^ = 
and ^ - - A w oscillatory phenomenon (Gibbs phenomenon). But in this case a closer solution 
for the given boundary condition is obtained. 

In analytical calculations, it is a common assumption to include only the first spatial eigenvalue 
Al. This one is associated with the longest escape timescale being the most important one for 
modelling escape of particles out of the galaxy. For numerical purposes, the computing time gives 
an upper limit to the possible number of eigenvalues. Furthermore, the spatial solution as performed 
in this paper has to match the formal solution for a spherical geometry within the limit of small 
eUipticity (e — > 0). ISchhckeiser et alj ( Il987b found as the spatial solution 



T(R, u) = y c„,—-Ji [y,„ — exp 



( Kly- 



2,,2 \ 



(37) 
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2 3 
f- 1 



Fig. 5. Spatial solution Eq.(l32]i in addition with weighting factors a, (Eq.(l20li) as a function of ^. 
As an illustration, we put - 4. For the solution represented by the solid line we included 200 
eigenvalues, whereas in the dashed solution only 10 eigenvalues are taken into account. 



corrected 



Schlickeiser et al 



which alr eady shows affinity to ou r solution Eg. (1321). The y,„ denotes the zeros of Ji(x) (here we 



(119871) ) and Ri describes the edge of the galaxy. Generally we can 



write for both solutions: 



T{X, m) = ^ aiZiiaX) sxp(-Aju) 



(38) 



The asymptotic behaviour of the solution Eq.(l32li is with Z(cjX) = ^ ^^^Ji (c,^) 

Z(c,f ) > — cos I C,f - -TT 1 . (39) 

The limit c-,^ oo satisfies a spherical symmetry with a numerical excentricity equal to zero (cf. 
Eq.(|23j). Furthermore with this limit, the semi-major axis c,^ goes to infinity. Performing the same 
limit to Eq.dlTl) with Z(c,X) = R'^^'^ J i{ciK), we get 



CiR^ao 1/1 

Z{ciR) > — cos \ciR - -TT 



(40) 



which is equal to Eq.([39]l. As a consequence we showed that our solution functions embody the 
spherical geometry within the limit of numerical excentricity equal to z ero. It is also straightforw ard 
to show that the weighting factors have the same structure as those in 
that the general solution is correct. 



Schlickeiser et al, 



(119871) so 



4.3. Momentum solution 

For the formal momentum solution, it is necessary to have a closer look at the spatial solution. As 
noted above, the eigenfunction expansion Eq.(fT9]l, in addition to Eq.(l20t. is indeed the best way to 
solve the spatial transport Eq. (fTT]) . Inserting this expansion into the convolution Eq.(fTOli. we can 
write 

/(r,/,) = 2 A,(r)7?,(p) (41) 

i 

with 

/-*QO 

Ri(p)= I duM(p,u)e-^^". (42) 
Jo 

Consequently the momentum solution Ri{p) obeys the ordinary diff'erential equation 

OpRi(p) - A]g{p)Rip) = -Q(p), (43) 
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which can be seen after multiplying Eq.(fT5]l by T(r, u) as given in Eq.(fT9]l and integrating over the 
convolution variable u from to oo. Each spatial eigenvalue /l^ enters this equation in the form of 
an inverse catastrophic loss time for particle escape out of the galaxy. Therefore Eq.(|43]) is called 
the 'leaky box' equation in momentum space. 

As the result of the eigenfunction Eq.(fT9ll. we have to solve one ordinary differen- 
tial equation for each spatia l eigenvalue instead of the partial differential Eq.lfTSl). Following 
Lerche & Schlickeiseii ( 119881) . we introduce some simplifying assumptions in order to find ana- 
lytic solutions of Eq.(l43ll. We take in Eq.dTSb 64 - 1, p{p) - p, and K(p) - (p/piY where pi is a 
normalisation value. These assumptions imply that the most important continuous loss process in 
elliptical galaxies at energies > 10 GeV is adiabatic energy loss due to a high galactic wind gradi- 
ent, whereas pion production losses are neglected because of the low number density of HI and HII. 
The fragmentation lifetime (oc 64) is independent of momentum, and the momentum dependence 
of the diffusion coefficient is defined by the parameter s - 2 - q, where q is the spectral index of 
the magnetic turbulence power spectrum. Therefore we get 



'1 



dp 

aAP\ 



mp;p — 



■ a^ip^Ri 



+ A1 



Ri = -Q(p)p-'PI 



(44) 



This ordinary difl'erential equation can be solved by a standard technique taking the finiteness 
of Rj at p ^ and p 00 into account. The formal mathematical solution is taken from 



Lerche & Schhckeiser 



19881). However, we found after correction of some minor mistakes: 



Riip) = (a2pir 



a+3 fi-[(3+a)/216i 
2s s{0\^,p,f'^- 



xp exp 



ySi +(/32+4iA,)'/2 

^ ' 

2 s 



xW 



a + 3 ^1 - [(3 -H fl)/2]j0i a + 3 OS? -h 4^,)i/2 
+ ■ . , . , P 



2s s(J3- 

rp 

I '^Po/?(r'2(Po)exp 
Jo 



4(A/)'/2 ' s ' 
/3i-082+4^.)i/2 



2s 



-Po 



xM 



+M 



a + 3 _^^i-[(3 + a)/2]/3i a + 3 (y6f-H4iA,) 



1/2 



2s 



s(J3\+Ail,iyi^ 



-Pi) 



>l/2 



a + 3 ^1 - [(3 -H fl)/2]j6i a + 3 (yet+^iA/) 

+ ^ : TT^, , 



2s i(y6f -h4iA,)1/2 



: I dpQ 
Jp 



/'(r'fi(Po)exp 



A -062+4iA,) 



1/2 



2s 



-Po 



xU 



a + 3 ^^i-[(3+ a)/2]/3i a + 3 (jSj + 4(A,) 



1/2 



2s 



s{0\+Aijfi)'l^ 



-Po 



(45) 



Here, we took for brevity 



Px = ail(a2p\), 

a = a\la2, 
^1 = a4l{a2p\). 



(46) 
(47) 
(48) 
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(49) 



The functions U and M denote confluent hypergeometric functions of first (Kummer) and second 
(Whittaker) order, respectively. 

For the special case of a (5-function injection of cosmic ray particles, 



Qip) = 6{p - Pinj), 

the solution Eq.(|45]l can be evaluated analytically. Therefore we get for p < 

r 



RAp) = {a2p\)-'- 



fl+3 , fi-[(3+a)/218i 
2s sOS2+^,.)i/2 



xp exp 



^ i 

2s 



^1/2 



fl + 3 - [(3 + fl)/2]/?i fl + 3 08f + 4(/r,)' 

+ ' ■ ■ . , , p- 



xpfnl' exp 



/3i +062+4^.y/2 ^ 



2,s 



inj 



fl + 3 ^1 - [(3 + fl)/2]ySi fl + 3 ifii+^ilJiY" , 
+ , , Pi„j 
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and for p > pi„j 



fl+3 , fi-[(3+fl)/2]/ii 

2.S i(j8j+l/(;)"- 



X;? exp 
xif/ 



2s 



fl + 3 ^1 - [(3 + fl)/2]^i fl + 3 C6^+4i/r,)'/^ ^ 
2s s(y82 +4iA,)i/2 ' s ' s 



xpfnl' exp 



/3i + (J3] + 4ifii) 



1/2 



2s 



-inj 



xM 



fl + 3 ^1 - [(3 + fl)/2]/?i fl + 3 08f +4iA,)'/2 
+ -:— , , p- 



2s sOe^ +4^,.)i/2 s 



inj 



(50) 



(51) 



(52) 



Because of the given structure of the solution Eq.(l45Tl. it is convenient to introduce the momen- 
tum value Pi, via 



s 



(53) 



This parameter indicates the momentum value above which the cosmic ray spectrum cuts off' ex- 
ponentially due to adiabatic losses (Pi) and escape losses (i/',). Now we can give some asymptotic 
behaviour for the solution: 



In the case of low values of the momentum {p « /?*), we find, according to lAbramowitz & Stegun 

mm . 



Up « p*) - p^'-\ 



(54) 
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Fig. 6. Cosmic ray spectrum (grey line) for a delta shape injection at p - p^j obtained by super- 
posing the contributions from individual modes (/ = 1,2, 3). The first mode Ri(p) modelling the 
longest escape losstime is responsible for the high energy cutoff whereas the slope of the power 
law is determined by the most efl'ective energy loss process (adiabatic loss). 



which is a flat power-law spectrum at very small momenta. For large arguments (p '» /?*) of the 
Kummer function in Eq.(|52l), we get 



exp 



2s 



1/2 



(55) 



In this situation, that the exponential cutoff" is not negligible. In the special case of dominating 
adiabatic losses, we take [f: » Aijji to get 



Ri[p » (i/y6i)'^') = exp 



s 



(56) 



Also in this case the exponential cutoff" dominates the spectral behavior. To summarize these re- 
sults we expect a ffat power-law spectrum at very low momenta, whereas the exponential cutoff' 
dominates at very high momenta. 

As an illustrative example a spectrum is modelled such that the result is a power law-spectrum 
with an exponent like the one observed from high energetic cosmic rays on Earth. This is demon- 
strated in Fig.®. Here particles are injected at a momentum value of p - pinj with a delta-shape 
injection function. Furthermore we assume an isotropic Kolmogorov turbulence model, i.e., the 
power law index of the turbulence is, q - 5/3. We found a steeper power-law spectrum as the 
predicted one from Eq.(l54l) over about two decades in momentum. In this regime the exponential 
cutoff" already plays a nonvanishing role. As the result of the small momentum dependence of the 
exponential term in Eq.dSTI) {s = 1/3 for Kolmogorov-Uke turbulence), a power-law spectrum over 
just two decades in momentum is obtained. Performing some delta-shape injections at increasing 
momenta with adequate normalisation values a power-law spectrum over a wide range in momen- 
tum is obviously possible. 

The parameter that aff'ects the final power-law index of the cosmic ray spectrum is the parameter 
a that describes the ratio of energy gains from the Fermi I process and the Fermi II momentum 
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diffusion process. To match the -4.8 spectrum from observations, the value has been chosen as 
a = 250. The adiabatic losses are assumed to be as strong as the Fermi II energy gains. In this 
context we need to remember that the solution function R/ip) has to be multiplied by Anp^ to 
obtain the number of particles at the momentum p: 

Niip) = 47Tp%(p). (57) 

The spectrum shown in Fig.© comprises the first 10 eigenfunctions, and the cumulative solution 
(grey line) is given by RuAp) - R\{p) + Riip) + •■• + Rxoip)- It can be seen that the spectrum 
is almost completely dominated by the first eigenfunction. Especially the cutofi" is dominated by 



Lerche & Schlickeiser 



this fir st Eigenfunction. That result has to be compared with, e.g.. Fig. 2 of 
(Il988h . where the situation in spiral galaxies is described. Both, elliptical and spiral galaxies have 
in common that the the cutoff" is described by the first Eigenfunction. The difi'erence between these 
two situations is that, at low energies, the higher Eigenfunctions describe the spectrum in spiral 
galaxies. The resulting power law is made by the sum of the single Eigenfunctions. In elliptical 
galaxies the situation is diff'erent, because the spectral shape is by oneself dominated by the first 
Eigenfunction. 

The reason for this different behavior compared to spiral galaxies is based on the loss time 
scales. While in the latter the escape loss processes dominate because of the small galactic height 
compared to the radial size, this is different in elliptical galaxies, where the smallest 'edge' of the 
confinement region is given by the semi-minor axis being larger than the galactic height in spirals. 
Therefore the dominating loss process in elliptical galaxies is adiabatic loss. 



5. Connection to observations 

The results of the solutions of the transport equations have been explained in the previous section. 
For astrophysical scenarios it is important to compare these with observations. Here we concentrate 
on the gamma-ray spectrum due to neutral pion decay and on the escaping cosmic ray spectrum 
from elhptical galaxies. 



5.1. Gamma-ray flux from pion decay 

The found solutions allow us to calculate the gamma-ray flux from pion decay. The pions are mainly 
produced due to proton-proton interactions, so we concentrate on the reaction p+p ^ p+p+xxn^. 
The neutral pion decays after the mean Ufetime T = 9xl0'^sin two gamma-photons. Other in- 
teraction channels can be treated analogously. Then the derivation is straightforward, but tedious. 
We have to convolve the high-energy proton spectrum, Eq.(l26b multiplied with Eq.(l52b. 
with the pion power of a single proton. The latter one is approximately given by (cf. 



Mannheim & Schlickeise 



]| 1 994I lschhckeiseJ[2002h : 



P.i7n, 7p, ^, 7j, 4,) = cy,T(^, 77, 4>)Ecj^°p{yp) x 6(y, - y'J'W [jp - y.hr] . (58) 

In this equation, Scr^^(yp) gives the multiplicity E for neutral pion production in p-p interac- 
tions and cr'^Ayp) the total cross-section. Within the high energy limit, this factor is given by (cf. 



Schhckeiseril2002): 



E 



s<(rp) = <pp(jp - (59) 
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The factor ythr, together with the Heaviside function, gives the threshold value for the proton energy 
to produce neutral pions (ythr'WpC^ - 1.22GeV). Therefore we get, for the pion source function. 



■„m„c^ J I 



xU 



dy„T(^,i],cf>)Ci,iy''/^exp 



2s 



A 1/2 



-(mpcYy 



fl + 3 ^1 - [(3 + fl)/2]ySi a + 3 + 4iA,)'^' 

+ — ■ — -7-z — , , (mpC) Y 



xcy„T{^, 77, 0, M)<^/ry, - 1)" '' X 5(y„ - y'J*)H [y„ - y.hr] 
The power-law factor p"^^ is due to Eq.dSTll and we used 



p 



(60) 



Cinj = 47T(mpCy^^ X (a2pl) 



n+3 , fi-[(3+fl)/2]/3i 



2s 



xPinV exp 
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xM 



a + 3 ^1 - [(3 + a)/2]^i a + 3 (j6^ + 4i/r,)'/2 



(61) 



Performing the integration we get 



4r2(^,/,,0)(rS, 



xCi„jrr'-^(r:/^ -!)«•« exp 



A +(/^^ + 4.A,)'^' 
2i 



Xf/ 



a + 3 , ^1 -[(3+fl)/2]^i fl + 3 (/??+4«A,)'^2^__ 4,^3 



2s s(j82 +4^,.) 1/2 



(mpc)- y 

(mpc)'y^ 



(62) 



It is obvious that the pion source function is proportional to the squared spatial distribution of the 
high-energy protons since pion production is a two-body process. Therefore the brightest luminos- 
ity is strongly correlated with the highest proton density in spatial coordinates because the mean 
lifetime of the neutral pion is extremely short so that the photons are produced very close to the 
p-p interaction region. 

With this result, the differential gamma-ray source function is given by (cf.li 



Schlickeiseril2002h 



/-•CO 

Qy(Ey,^,r],^) ^2 I dy„- , 



(63) 



where ymin - + '^j§- denotes the minimum energy of the gamma-ray photons after the decay 
of the neutral pion. We can state that the differential gamma-ray source function in the high-energy 
regime {Ey » m„oc^), where the spectrum shows a power-law dependence, is approximately given 
by 



xO.53 , 



(OT,,c)Vp' 



xU 



fl + 3 , ^1 - [(3 + fl)/2]^i fl + 3 4ipd"^ 4,/3 



-(mpcYyl'' 



(64) 



2s s(/3]+4iJ/iy/2 s s 
but the integral Eq.(l63Tl cannot be performed analytically. 

Consequently, as discussed before, we assume a power-law spectrum over a wide range in momen- 
tum until a maximum value of 7p, max'W/jC^: 

N(yp) = Noy-'H [yp, - yp] . (65) 
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After performing the convolution (cf. Eq.(|60ll) the pion source function reads 

2 ^ 

QAJ., ^, V,4>)-\ Li!±^^^^-(4.-l)/3(^4/3 _ 1)0.53^ _ 3/4J ^ [^3/4 _ ^^^^ 

i nijic ^ J L J 

From this calculation we get due to Eq.(l63]l in the high-energy limit ( yjyl - 1 ^ y„ and {yl^^ - 
1)0-53 ^ y^^"'^"'^"') for gamma-ray energies above 130 GeV 



(67) 



as the gamma-ray source function. Thus the differential gamma-ray source function at high energies 
(Ey » m^c^) is given by the relation 

QyAEy, ^, ^, 4') ^ E--'' = (68) 

after Eq.(l66b. Taking an injection power-law index z = 2.8, we find that Qy „o(Ey, ^, rj, (p) cc Ey^-^^, 
which is slightly harder than the injected proton spectrum. Therefore an observed gamma-ray flux 
from an elliptical galaxy gives constraints to the high-energy cosmic ray spectrum within that 
object since neutral pions are produced in the whole confinement volume. Therefore we expect an 
elliptical galaxy be an extended source of gamma rays. 



5.2. Escape-spectrum of cosmic rays 

We next derive the spectrum leaking out of elliptical galaxies. Therefore the mean free path Amfpip) 
for single scatterin g events between charged particles and plasma waves is needed. According to 



Schhckeisen (l2002h . this characteristic length is given by 

^mfp(p) = H^Il^I ^ Ao.mfpp^-'^, (69) 



where K(r, p) and v denote the spatial difiiision coefficient and the particle velocity respectively. 
Again the parameter q gives the momentum dependence of the turbulence spectrum. Here we use 
V - c for highly relativistic cosmic rays. FiguredTJi gives the chosen geometry for calculating of the 
escape spectrum. We assume that cosmic rays will escape out of the galaxy, if their mean free path 
is larger than the distance to the edge of the galaxy. This border value is given by (cf. Ch. (l4.3b ) 

r,(^, 77, 0) = ^(^77)2 + (^2 _!)(!_ ^2), (70) 

All particles in the gray-shaded area will escape the galaxy, if rg(^, 77, (p) - d\ < A^f^ip), where we 
neglect a geometrical factor of order unity. The width of this spherical shell is constant in every 



direction so that fe -\/^2 _ 1 - d\ = fe^ - d2, as indicated in Fig.©. The rate of the escaping high 

energetic particles is approximatively given by the inverse escape time 

K{r, p) ^ cA^fp(p) 

rH^,T],4>) 3/2[(^77)2 + (^2-l)(l-772)]' 

where we neglect the geometrical factor for the direction of the trajectory of cosmic ray particles. 

Therefore the total number of particles escaping out from the gray-shaded area into the intergalactic 

medium per unit time and per unit momentum is given by 



C/li 

3/2 [(^77)2 + (^2-l)(l- 772)]^ 



X// [A^f,(p) - (r,i^, 77, 0) - d,)] ^^^^^^^-^^fdfdr^dcp. (72) 
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Fig. 7. Schematic view of the geometry for calculating the escape-spectrum from elliptical galaxies. 
Particles within the gray-shaded area leave the galaxy by chance, if their mean free path A^fpip) is 
longer than the distance to the edge of the galaxy. 



The upper boundary in the /-integral gives the radius of the galaxy with respect to the vari- 
ables ^ and 77. As an illustration we choose a constant cosmic-ray distribution like Eg. ( II 06b with 
Eq. dlOTK see Appendix B). This one is independent of the variables 77 and 0. More realistic particle 
distribution functions can be treated analogously. Furthermore, we take a constant value of ^ for a 
dedicated elliptical galaxy. The Heaviside function H [/imfpC/') - '?> 0) ~ "^i)] impUes that only 
high energetic particles with /Imfp(p) > fgi^, 77, 0) - di can leave the galaxy. 
Under these conditions we can convert Eq.(l72]l into 

(73) 

The general solution of this integral is given by 



16n^cp^Al ip) cot-i( Vi^) V- 



(74) 



From this calculation we recognise that the escaping spectrum we observe depends on the shape of 
the galaxy given by the choice of the parameter ^=constant. Inserting Eq.(l69t into Eq.(l74l). we get 



16 2 2 I 
oc — TT C/ig ,^fp 2^ OiTii^ = const)Ri(p)p 



6-2q 



(75) 



Under the assumption that Rges(p) is proportional to /? we find for a Kolmogorov turbulence 
model (q - 5/3) a power-law index of -2.13 for the escaping spectrum. This spectral behaviour is 
harder than the intrinsic spectrum, which has a spectral index of -2.8 (cf. Ch.( l4.3b ). 

Such an effect can be easily understood by the high energetic particles leaving the galactical 
confinement region more frequently than low energetic ones. As a consequence there are two pos- 
sibilities explaining the observed high energy cosmic ray spectrum above the 'knee' with spectral 
index of -3.1 assuming that elliptical galaxies provide a significant amount to the overall high- 
energy cosmic ray flux. First, the intrinsic spectrum in elliptical galaxies may be steeper then that 
one we have chosen here. In this context the measurement of the gamma-ray flux reaching the Earth 
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from such an object would give us interesting constraints. Second, if the intrinsic spectrum in ellip- 
tical galaxies has nearly the same dependence as the one we measure here on Earth {N(p) oc p^^-^), 
the transport of cosmic rays after escape from elliptical galaxies depends on energy. In these two 
scenarios it is possible to explain the high-energy component with our model. 

Note that, due to the previously given arguments, the highest energetic particles cannot be 
confined within the galactical volume. In our examination the maximum energy of confined charged 
particles within elliptical galaxies is given by the relation 

fe^ = y^mfpip) = \mfpP^'''- (76) 

For reas onable values of g iant elliptical galaxies (/lmfp(p = 10 GeV c"' ) = 10''' cm, fe^ - 5 x 10^^ 
cm, (cf. Schlickeisetll2002 )) we get from Eq.(|76l) a value of about 10^' eV for cosmic ray protons. 



Above this energy value, it is generally impossible that elliptical galaxies confine cosmic rays 
within their volume. 



6. Conclusion 

We showed that an analytical treatment of cosmic ray transport in elliptical galaxies based on 
the diffusion approximation is possible in general. The formal solution we found, combined with 
appropriate boundary conditions and source functions, can be used to study transport processes in 
elliptical galaxies. This model is valid for the complete physical parameter space as long as the 
diffusion approximation holds. 

The first test case where we applied our model to is a jet-like injection shape like in M87. 
As we can separate our problem into spatial and momentum problems, this test case probes the 
spatial problem. For this model we found that for very long timescales cosmic rays are distributed 
throughout the galaxy almost isotropically. 

We also provided a test case for the momentum problem. Under the assumption of a resulting 
power-law spectrum matching the observed power law-spectrum on Earth, we could identify the 
governing eigenfunctions and therefore the governing processes for cosmic ray acceleration and 
momentum diffusion in elliptical galaxies. It turned out that in elliptical galaxies adiabatic losses 
are responsible for the high energy cutoff, whereas the slope of the spectrum is given by the ratio 
of the strength of Fermi I and Fermi II processes. As a result of the basic possibility of explaining 
this power law-with physical parameters that might be found in M87, our calculation gives rise to 
the theory of M87 as a source of ultra-high energy cosmic rays. In this context the gamma-ray ffux 
from pion decay and the escaping cosmic ray spectrum are essential for a more adequate view of 
cosmic ray astrophysics. We find that the gamma-ray spectrum with a power law index of -2.69 
is a bit harder than the intrinsic cosmic ray spectrum (power law index -2.8). In the context of 
understanding the sources of high-energy cosmic rays, the ffux of charged particles escaping from 
the confinement region of elliptical galaxies is also very interesting. For the same intrinsic spectrum 
as above, we find -2.13 as the spectral index for escaping particles. Again the spectrum is harder 
than the value within elliptical galaxies. 

Our model may also serve as a testbed for some more advanced questions: 

- Are elliptical galaxies, and especially AGN within these objects, reasonable sources of the 
ultrahigh energy cosmic rays? 
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- What is the total flux leaving a single elliptical galaxy with respect to a measured source of 
cosmic rays in elliptical galaxies? 

- What is the total cosmic ray flux we can expect to measure on Earth in the special case of the 
nearby active elliptical galaxy M87 or Centaurus A? 

Clearly, to answer these questions, a detailed quantitative analysis of the properties of cosmic ray 
transport within elliptical galaxies is required. Therefore we need to determine from observations 
the exact physical parameters of the interstellar medium and of the emission processes within the 
source. With our analytical model at hand, the task of constraining the possible parameter space 
gets easier. Especially by assuming that UHECR are coming mostly from sources like M87, we 
may be able to derive the total cosmic ray content of this type of elliptical galaxies. 

Nevertheless this work is considered as a starting point for more sophisticated (numerical) 
models. Another interesting point is the expansion of this model to cosmic ray electrons. While the 
spatial description may be derived analogously, we need to add the physical processes involved in 
the transport of electrons in momen tum space. Synchrotron losses and the inverse Compton effect 



will then play a major role (see, e.g.. lCasadei & Bindil2004l) . The output of this model can be easily 



tested with the radio data of elliptical galaxies since the electrons provide the main contribution to 
this radiation. 
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Appendix A: General solution of Eq.(|25t 

The general solution of Eq.dZSll can be obtained by consecutive separation of variables. First we 



use 

r(^,77,,^,M)=X(^,77,,^)f/(M), 

leading to the following dififerential equations with -k^ as the separation constant: 

au 



d_ 



dr] 



The factor c is given by 



orj 

+ c^if - rf)X = 0. 



Equation ( |78] i is easily solved by 



U{u) - exp(-^ m), 



(77) 



(78) 



(79) 



(80) 



(81) 



where we put the integration constant to equality without loss of generality. This result serves 
an a consistency check to the eigenfunction expansion Eq.(fT9]l. The separation constant k will be 
calculated by use of spatial boundary conditions. 
Eq.(|79]l can be solved through another separation of variables. With 



■ A{^,ri)B{<t,), 



(82) 
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one gets 



1 d'^B _ (^^- 1)(1 -77^) 
"S502 ~ (^2 - 772) X A 



d_ 

drj 



OT] 



The separation constant is set to be +m . Therefore we get 

—5- + rn^B = 0, 



and 



ie - i)(i d 



(f-rf-)xA 



(r - 1) 



dA 
5? 



d_ 
drj 



(1-1^) 



+ [(f--l)(l-T]^)c^-m^]A^Q. 
The solution of Eq.([84]i is given by 

= sin(m0) + cos(m0), 
where we again used integrating constants equal to 1. The last separation 

A(^,77)=7?©5(77) 

solves Eq.(|85]l. which can be written by 



d_ 



d_ 
di] 



1 



I— 

1 



[f-l 1-772 
With the separation constant +/?,„„ we get 



A =0. 



and 



d_ 

d_ 

drj 



07] 



e-i 



+ {l3,n„-r]'c'-j—^]S =0. 



(83) 



(84) 



(85) 



(86) 



(87) 



(88) 



(89) 



(90) 



It is notable that the ^ and 77 solution functions satisfy similar differential equations. From 



Abramowitz & StegunI (Il972h we find 

1 



^m«(^' ^) 



00' 
r=0,l 



(2m + r)! 



X 



m/2 



(91) 



as the solution for the ^ dependence of Eq.dZSTl. It is zl''\c^) = -^^■^«+i('^?) for the first kind 
(7:1 = 1) and zll'\c^) - .J^Y^^i(c^) for the second kind (p - 2). J and F are Bessel functions 
of the first and second kind respectively. For our model we can neglect the solution of second 
kind because of physical motivated spatial boundary conditions (see Appendix B). It is clear that 
n - m + r. So in Eq.(|9TTl is equal to 1 . The mark ' connected with the sum means a summation 
either over even or over odd values of r. The latter ones as well as the weighting factors d"'" have 
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also been adjusted to spatial boundary conditions. 
Eq.(|90]l has the formal mathematical solution given by 

oo' 

5S(c,/7)= ^^rWi^^+.W, (92) 

r=0,I 

and 

co' 

5£i(c,/7)= 2^r(c)e;;;,.w- (93) 

r=0,l 

Here Pm+r(v) and denote associated Legendre polynomials of first and second kind respec- 

tively and d'^"(c) are weighting factors. Again ' means a summation over even or odd values of r. 
We can neglect the solution function of second kind with respect to spatial boundary conditions 
(see the consistency checks). 

Under these conditions we can write the formal mathematical solution as 
m,7],cf>,u) = YjR^^l(c,OxSlllic,i])xcos(mcp)xexp{-k\)^ 

m,n 

m,n 1 Lr=0, 1 

^ 7 (2m + r)! 



r=0,l 



x<"(c)P;;;+,(77) x cos(m0) X exp(-k^u)] . (94) 

Appendix B: Effect of spatial boundary conditions to tlie general solution Eq.((26) 

As an illustrative example we calculate the solution depending on the following boundary condi- 
tions. First we define one for the variable ^. Remembering, that is equal to the semi-major axis 
a we assume 'free escape' boundary conditions in the form 

r(^ = ^„77,0,M) = O. (95) 

Hence the particle distribution function T{^, rj, <p, u) is set to be zero at the maximum size f^c of 
the confinement region of cosmic ray particles being a convenient way to model leakage out of the 
galaxy. This region is assumed to be larger than the (optical) size of the galaxy (/^g). Furthermore 
we assume periodical boundary conditions for 77 and ^: 

T(^,j]=l,^,u)^T{^,T] = +l,(f>,u) (96) 

Ti^, = In, u) = T{^, T],(p = 0, u) (97) 

Finally we take 

r(^,77,0,M = O) = ^o0(^g-^) (98) 

for a given distance of the foci (2/) as a simple model to explain a constant source function over the 
whole size of the galaxy. Here we neglect the (p- dependence (taking a 0- dependence into account 
is straightforward). As a consequence we set m = in Eq.(l26ll leading to n = r (cf. Appendix A). 

In order to match Eq.(|96ll only even solution functions are practical. This is done by P^i.{jf). 
Therefore we neglect Q%rj) having singularities at the points rj - \ and 77 = -1. From Eq.(|95l) we 
recognise that J^.^\_{c^) - Oif ^ - ^c- Let be the zeros of J^.^\_{^), then we find 

C = — , kri^ ■ (99) 

5c fcj 
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As a result we write 

- - ■ 

This means that we have to sum over all zeros y„ aditionally. The last boundary condition Eq 
is useful to determine the weighting factors c/,-, and dri respectively. From 



(100) 



2 /■=(), 1 ^ri 



2 dnP°(v) 



r=0,l 



(101) 



We see that we can cancel the factors c/,,. 

Adjusting Eq.(l26]l to the boundary condition Eq.(l96b it becomes clear that only even values of 
r match the required periodicity. Using the orthonormality relation for Legendre polynomials 



X 



' 2 



2r+ 1 



the RHS of Eq.lfTOTTi evolves to be 



(ri)dTj. 



(102) 



(103) 



Performing the integral it is obvious that only r —Q gives a non-vanishing value for j^'^ F^irfjclri - 
2. Consequently we set r = r = 0. As a result all summations over r disappear so there is only 



left. Here we use the orthonormality relation for Bessel functions, i.e. 

.2 



(104) 



(105) 



to calculate the weighting factors. The Bessel function of second kind diverges at ^ = and is 
therefore not appropriate for the general mathematical solution (cf. Appendix A). Finally we find 
(we substitute a/ = J^r ^ di) 



£,(r) = T(^, «) = ^ "' "^"^5 [^'j') 



with the weighting factors 



(106) 



(107) 



f h(3'.)l 

Again we see that the general solution of Eq.dZSll with respect to spatial boundary conditions is 
indeed of the form of the eigenfunction expansion Eq.(fT9]l with Eq.(|20]i. 
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